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Abstract 

Spatial extent — the possibility that what happens at point A cannot immediately affect what 
happens at point B — is a complicating factor in mathematical biology, as it creates the opportunity 
for spatial non-uniformity. This non-uniformity must change our understanding of evolutionary 
dynamics, as the same organism in different places can have different expected evolutionary out- 
comes. Since organism origins and fates are both determined locally, we must consider heterogene- 
ity explicitly to determine its effects. We use simulations of spatially extended host-pathogen and 
predator-prey ecosystems to reveal the limitations of standard mathematical treatments of spatial 
heterogeneity. Our model ecosystem generates heterogeneity dynamically; an adaptive network 
of hosts on which pathogens are transmitted arises as an emergent phenomenon. The structure 
and dynamics of this network differ in significant ways from those of related models studied in the 
adaptive-network field. We use a new technique, organism swapping, to test the efficacy of both 
simple approximations and more elaborate moment-closure methods, and a new measure to reveal 
the timescale dependence of invasive-strain behavior. Our results demonstrate the failure not only 
of the most straightforward ("mean field") approximation, which smooths over heterogeneity en- 
tirely, but also of the standard correction ("pair approximation") to the mean field treatment. In 
spatial contexts, invasive pathogen varieties can prosper initially but perish in the medium term, 
implying that the concepts of reproductive fitness and the Evolutionary Stable Strategy have to 
be modified for such systems. 



1 



I. INTRODUCTION 



Mathematical modeling of biological systems involves a tradeoff between detail and 
tractability. Here, we consider evolutionary ecological systems with spatial extent — a compli- 
cating factor, as it implies the possibility of heterogeneity. Analytical treatments of spatial 
systems typically treat as equivalent all configurations with the same overall population 
density, the same allele frequencies, the same pairwise contact probabilities or the like. Gen- 
erally, one seeks a simplified analytical model, which coarse- grains "microstates" (complete 
specification of each organism) to "macrostates" (characterized by quantities like average 
densities), allowing one to make useful predictions about the model's behavior [1, 2]. Cor- 
rections to simple coarse-grainings can quickly generate an overbearing quantity of algebra. 
Before exhausting ourselves with ever-more-elaborate refinements, it would be useful to have 
some understanding of when a particular series of approximations is doomed to inadequacy. 

In this article, we review the context in which commonly-used coarse-grainings can be 
expected to fail at capturing the evolutionary dynamics of an ecosystem, and in addition 
we provide a novel, direct demonstration of that failure. The fundamental issue is spatial 
heterogeneity, a long-recognized concern for mathematical biology [3, 4]. When does spatial 
heterogeneity significantly impact the choice of appropriate mathematical treatment, and 
when does a chosen mathematical formalism not capture the full implications of spatial vari- 
ability? We show that one can test a treatment of heterogeneity by transplanting organisms 
within a simulated ecosystem in such a way that, were the treatment valid, the behavior of 
the ecosystem over time would remain essentially unchanged. We demonstrate situations 
where the behavior changes dramatically and cannot be captured by a conventional treat- 
ment. The complications we explore imply that short-term descriptions of what is happening 
in an evolutionary ecological model can be insufficient and, in fact, misleading, with regard 
not just to quantitative details but also to qualitative characteristics of ecological dynamics. 

Many modeling approaches in mathematical biology which appear distinct at first glance 
turn out to be describing the same phenomenon with different equations [5, 6], e.g., by a 
different choice of coordinate system [7]. What matters for our purposes is not so much 
which technique is chosen, but whether the underlying assumptions do, in fact, apply. 

It will be convenient to introduce a small amount of jargon: "mean-field theory" is a 
term from statistical physics [8, 9] which has been adopted in ecology [10-12], referring 
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to an approximation in which each component of a system is modeled as experiencing the 
same environment as any other. This imphes that the probabihty distribution over all 
possible states of the system factors into a product of probability distributions for individual 
components. In a mean-field treatment of a magnet, for example, each atom experiences 
the same magnetic field. An analogue in population genetics is the assumption that a 
population is panmictic. That is, if a new individual in one generation has an equal chance of 
receiving an allele from any individual in the previous generation, then we can approximate 
the ecosystem dynamics using only the proportion of that allele, rather than some more 
complicated representation of the population's genetic makeup. Modeling evolution of that 
population as "change in allele frequencies over time" (per, e.g., [5, 13]) is, implicitly, a mean- 
field approximation [14]. The mean- field approximation is also in force if one postulates that 
an individual organism interacts with some subset, chosen at random, of the total population, 
even if the form and effect of interactions within that subset are complicated (as in, e.g., 
[15, 16]). 

It is well known that real species are not necessarily panmictic. However, many treat- 
ments which acknowledge this are still mean-field models. The textbook way of incorporating 
geographical distance into a population-genetic model is to divide the system into local 
subpopulations, "islands," connected via migration [17-19]. Within each subpopulation, 
distance is treated as negligible, and organisms are well mixed [4, 20]. This approach as- 
sumes a single distance scale below which panmixia prevails [21], and it relies on well-defined 
boundaries which persist over time [20]. Furthermore, the connections among subpopula- 
tions are frequently taken to have the topology of a complete graph, i.e., an organism in 
one subpopulation can migrate to any other with equal ease [4, 18-20]. In this case, each of 
the N subpopulations do experience the same environment, to within one part in A^. Thus, 
the mean-field approximation is in force at the island level, and the island model incorpo- 
rates spatial extent without incorporating full spatial heterogeneity. For real ecosystems 
[21-23], or in sufficiently interesting models, one or more of these simplifying assumptions 
can fail. More complicated population structures require more sophisticated mathematical 
treatments of evolution, a fact which has real-world implications for practial issues like the 
evolution of drug-resistant diseases [24]. 

We make these issues concrete by focusing on a specific model of ecological and evolu- 
tionary interest. For definiteness, we take a model of hosts and consumers interacting on a 
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2D spatial lattice. Each lattice site can be empty (0), occupied by a host (H) or occupied 
by a consumer (C). We use the term consumer as a general label to encompass parasites, 
pathogens and predators. Hosts reproduce into adjacent empty sites with some probability g 
per site, taken as a constant for all hosts. Consumers reproduce into adjacent sites occupied 
by hosts, with probability r per host; sometimes r is fixed for all consumers, but we shall 
also consider cases in which it is a mutable parameter passed from parent to offspring. Hosts 
do not die of natural causes, while consumers perish with probability v per unit time (leav- 
ing empty sites behind). Because consumers can only reproduce into sites where hosts live, 
the effective graph topology of reproductively available sites experienced by the consumers 
is constantly changing due to their very presence — it is an adaptive network [25-27]. This 
model is a fundamentally nonequilibrium one: an approach to a steady state is not at all 
guaranteed, and even when a quasi-steady-state behavior emerges, as we shall see, it is a 
consequence of fluctuations over space and time. 




FIG. 1: Snapshots of a simulated host-consumer ecosystem on a 250 x 250 lattice, taken at intervals 
of 100 generations. Consumers are dark gray (red online), hosts are light gray (green online) and 
empty space is left white. The simulation began with a single consmner at the center of the lattice, 
which gave rise to an expanding front of consumers. The first image in this sequence shows the 
state of the ecosystem 100 generations into the simulation. Hosts which survive the consumer wave 
recolonize the empty sites, leading to pattern formation. Here, the host growth rate is g = 0.1, the 
consumer death rate is i; = 0.2 and the consumer transmissibility is fixed at r = 0.33. 

Several different types of biological interactions can be treated by this modeling frame- 
work. Hosts could represent regions inhabited by autotrophs alone, while consumers repre- 
sent regions containing a mixture of autotrophs and the heterotrophs which predate upon 
them, for example [28]. Alternatively, host agents could represent healthy organisms, while 
consumers represent organisms infected with a parasite or pathogen. Thus, host-consumer 
models are closely related to Susceptible-Infected-Recovered (SIR) models, which are epi- 
demiological models used to understand the spread of a disease through a population. SIR 
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models describe scenarios in which each individual in a network is either susceptible (S) to 
a pathogen, infected (I) with it, or recovered (R) from it; susceptible nodes can catch the 
disease from infected neighbors, becoming infected themselves, while nodes which have be- 
come infected can recover from the disease and are then resistant against further infection. 
Susceptible, infected and recovered individuals roughly correspond to hosts, consumers, and 
empty cells, respectively. An important difference between host-consumer models and epi- 
demiological models concerns the issue of reinfection. To wit: in the host-consumer model, 
an empty site left behind by a dead consumer can be reoccupied by another consumer, 
but only if a host reproduces into it first. Other research has considered models where 
R[ecovered] individuals can also become I[nfected], with a different (typically lower) prob- 
ability than S [usceptible] ones, thereby incorporating imperfect immunity into the model 
[29, 30]. The degree of immunity is independent of geography and the environment of the 
R[ecovered] individual, unlike reoccupation in the host-consumer model. Another appli- 
cation is illustrated by the Amazon molly, Poecilia formosa, which is a parthenogenetic 
species: P. formosa, all of which are female, reproduce asexually but require the presence 
of sperm to carry out egg development. (This kind of sperm-dependent parthenogenesis is 
also known as gynogenesis.) P. formosa are thus dependent on males of other species in the 
same genus — usually P. mexicana or P. latipinna — for reproduction. Because P. formosa 
do not incur the cost of sex, they can outcompete the species on which they rely, thereby 
possibly depleting the resource they require for survival, i.e., male fish [18, 31]. Thus, hosts 
could be regions containing sexual organisms, with consumers standing for areas containing 
both sexual and asexual individuals [18]. 

This host-consumer model displays waves of colonization, consumption and repopulation. 
Hosts reproduce into empty sites, and waves of consumers follow, creating new empty regions 
open for host colonization. Therefore, clusters of hosts arise dynamically [32-35], a type of 
pattern formation which can separate regions of the resources available to pathogens into 
patches without the need for such separation to be inserted manually. Figure 1 illustrates 
a typical example of this effect. This is a specific example of the general phenomenon of 
pattern in nonequilibrium systems [14]. Consumers are ecosystem engineers [36, 37] which 
shape their local environment: an excessively voracious lineage of consumers can deplete the 
available resources in its vicinity, causing that lineage to suffer a Malthusian catastrophe 
[18, 32, 38-41]. Because the ecology is spatially extended, this catastrophe is a local niche 
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annihilation, rather than a global collapse [42]. A ravenous mutant strain can successfully 
invade in the short term but suffer resource depletion in the medium term, meaning that 
averages taken over long numbers of generations yield a relatively low virulence [10, 43]. 
This implies that an empirical payoff matrix or reproduction ratio will exhibit nontrivial 
timescale dependence [10, 44, 45]. 

Stochastic departures from uniformity — fluctuations — of this kind have been found to be 
important in the context of phase transitions, shifts from one "phase" or regime of behavior 
to another, such as the change of a substance from liquid to gas. Indeed, there exist sig- 
nificant and useful connections between phase transitions in physics and the present model. 
In certain situations, the process of pathogen propagation through the host population dis- 
tributed in space can be mapped onto a percolation problem. A topic widely investigated 
in mathematics, percolation deals with movement though a matrix of randomly placed ob- 
stacles. A prototypical percolation problem is a fluid flowing downhill through a regular 
lattice of channels, with some of the lattice junction points blocked at random. The key 
parameter is the fraction of blocked junction points. If this fraction is larger than a certain 
threshold value, the fluid will be contained in a limited part of the system. However, if 
the blocking fraction is below the threshold, the fluid can percolate arbitrarily far from its 
starting point. This is a phase transition between a phase in which fluid flow can continue 
indefinitely and one in which flow always halts. Similar issues arise when a pathogen prop- 
agates by cross-infection through a set of spatially arranged hosts. Sufficiently many hosts 
in mutual contact are required for the pathogen to propagate successfully. Pathogen strains 
therefore survive or die out over time depending on whether percolation is or is not possible 
[12, 46-50]. 

The mathematical connection between pathogen-host and percolation problems can pro- 
vide insight into the difficulties in analytical treatment of the biological problem. Spatial 
heterogeneity gives rise to failure of traditional analytic treatments of percolation and a need 
for new methodologies. Since the pathogen problem maps onto the percolation problem un- 
der some circumstances, the same analytic problems must arise in the biological context. 

Quite generally, in systems which have spatial extent, stochastic variations close to a 
phase-transition point cause mean-field approximations to fail [30, 51]. Moreover, even 
standard types of corrections to mean-field approximations perform poorly or fail entirely 
[11, 52]. The question becomes whether models of interest contain such troublesome regimes. 
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and how likely we are to encounter them. 

We will show that numerical investigation of our host-consumer model indicates that the 
consumer extinction transition, when the transmissibility r becomes just large enough that 
the consumer population sustains itself, lies in the directed percolation [30, 53-56] universality 
class. A similar result has been reported for related models [46, 47]. The directed percolation 
universality class is a large set of models, all of which exhibit a phase transition between two 
regimes of behavior, and all of which behave in essentially the same way near their respective 
transition points. The scenario of fluid flow through a random medium considered above 
is a classic example of a directed percolation-class model, but many others exist as well 
[30, 53]. The critical exponents describe how properties of the modeled system vary over 
time or as a function of how far the control parameter is from the critical point. They are the 
same for all systems in the universality class. Other universality classes exist as well, with 
different classes having different quantitative values for the critical exponents. Identifying 
to which universality class a system belongs enables us to study a complicated phenomenon 
by examining a simpler representative of its class instead. 

The inadequacy of the mean-field approximation implies a need for improved treatments. 
To include heterogeneity, as a first correction, we can include the pairwise correlations 
between sites. The details of the resulting calculations depend on the particular system 
being studied, but all these specific implementations collectively fall under the heading of pair 
approximations. This scheme is a common way of building a mathematical treatment which 
acknowledges spatial heterogeneity, by incorporating the idea that the local environment of 
an organism can have a different character and composition than the ecosystem as a whole. 
This is only an approximation, because it does not reflect the existence of patches (empty 
areas, host areas, and pathogen areas) and other population structures which involve more 
than two sites. We will see in §IV that pair approximations fail to capture what happens 
when an invasive species is introduced to a host-consumer ecosystem. We will show two 
reasons for this failure, phase transitions and long-times cale effects. Furthermore, these 
properties imply that including triple, quadruple or even higher order correlations will also 
be insufficient. The failure is due not to the details of any particular calculation but rather to 
an assumption inherent in these approximations. Thus, there exists a need for new analytical 
tools in spatial evolutionary ecology. 

The failure of analytic approaches to address a particular problem often leads one to 
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consider alternative problems that can be addressed using the analytic approach. However, 
such a change, by its very nature, may eliminate effects that are important in the original 
problem, and consequently the solution does not give the correct insight into them. For 
example, there are host-pathogen systems defined on networks for which pair-approximation 
schemes can work. Prior modeling efforts have considered epidemics on adaptive networks, 
where the spread of the disease through the network changes the connections of the network 
[57-64]. In such models, if a susceptible node has an infected neighbor, it can break that 
connection by rewiring to another susceptible node. The key point is that the new neighbor 
is chosen at random from the eligible population. This choice of rewiring scheme is exactly 
what makes a pair approximation work, because it eliminates higher-order correlations in 
the system [58]. In our system, by contrast, hosts can form new connections by reproducing 
into empty sites, but these contacts can only connect geographically proximate individuals. 

While the presence of a nonequilibrium transition point indicates that traditional analysis 
techniques fail, it raises the possibility that new tools from the theory of phase transitions 
[30, 51] will become applicable. 

In §11, we start to address the issue of how mutation and natural selection operate together 
in a spatial host-consumer ecosystem. This leads us into the study of phase transitions and 
the relationship between ecological dynamics and percolation problems. In §111, we explore 
how to define fitness in the nonequilibrium context. In §IV, we will make direct tests of 
assumptions used in developing mathematical treatments for spatial evolutionary ecology. 
Imphcations for questions of wider scope are the subject of §V. 

II. SUSTAINABILITY AND PHASE TRANSITIONS 
A. Evolution of Transmissibility 

We investigate evolution in the spatial hostconsumer ecosystem through simulation and 
analytic discussion. If the transmissibility r is made a heritable trait, passed from a consumer 
to its offspring with some chance of mutation, what effect will natural selection have on the 
consumer population? Figure 2(A) shows the average, minimum and maximum values of 
the transmissibility r observed in a population over time. The average r tends to a quasi- 
steady-state value dependent on the host growth rate g and the consumer death rate v; if 
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the simulation is started with r set to below this value, the average r will increase, and 
likewise, the average r will decrease if the consumer population is initialized with r over 
the quasi-steady-state value. Even when the average r has achieved its quasi-steady-state 
value, the population displays a wide spread of transmissibilities whose extremes fluctuate 
over time [28]. 




FIG. 2: (A) Minimum, average and maximmn transmissibility r for a consumer population over 
time, with g = 0.1 and v = 0.2. (The mutation rate is /u = 0.255 and the step size is At = 0.005, 
as was used in reference [28].) The average r tends to a quasi-steady-state value dependent on 
g and v; if the simulation is started with r set to below this value, the average r will increase, 
and likewise, the average r will decrease if the consumer population is initialized with r over the 
quasi-steady-state value [28]. (B) Standard deviation of the transmissibility, CTt-j for a consumer 
population {g = 0.1, v = 0.2) measured at different mutation rates, from n = 0.003 to // = 1.0. 

In a well- mixed ecosystem, the average r of the population will tend to 1, maximizing the 
reproductive rate of the individual consumer. The observation of a quasi-steady-state value 
below 1 is an important result. This is the first breakdown of the mean-field approximation, 
and it indicates the inapplicability of traditional assumptions about fitness optimization, 
with implications for the origins of reproductive restraint, communication-based altruism 
and social behaviors in general [10, 28, 42, 44, 45, 65]. 

To gain insight into the upper and lower extremes seen in Figure 2(A), we study the 
simpler case of fixed r. We fill the lattice with hosts and inject a single consumer with a 
selected value of r; then, we observe how long the descendents of that consumer persist as a 
function of r. The consumer population does not persist when r is either too low or too high. 
Figure 3 shows the probability that a consumer strain will survive for a substantial length 
of time (2000 generations) after injection into a lattice filled with hosts. This probability is 
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hump-shaped, with an asymmetric plateau bounded above and below by cutoffs. 



B. Lower Threshold: Percolation 



1.0 




Transmissibility x 

FIG. 3: Probability of a consumer strain surviving 2000 generations after injection at a single point 
in a 250 x 250 lattice filled with hosts. (Computed with 1000 runs per point. Host reproduction 
rate g = 0.1, consumer death rate v = 0.2.) Vertical dotted line shows the sustainability threshold 
found through the mean-field approximation. 

The lower edge in Figure 3 is a transition between sustainable and unsustainable trans- 
mission of consumers through the lattice. We expect and will demonstrate that the system 
maps onto a directed percolation problem. We confirm that this is a directed percolation- 
class phase transition by finding the power-law growth of the population size, averaged over 
an ensemble of 10^ simulation runs. As shown in Figure 4(A), the exponent of this power 
law is the characteristic value for directed percolation, 0.230 [30]. If the host growth rate 
parameter is made negligible {g ^ 0), then the ecosystem behavior changes. The transition 
instead belongs to the dynamic percolation universality class, examples of which have been 
seen in some epidemic models [30, 47]. We can construct a prototypical example of the 
dynamic percolation class by returning to the fiuid-fiow problem considered earlier and, in- 
stead of requiring that the fluid only flow downhill, allowing it to flow isotropically through 
the lattice. 
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FIG. 4: (A) Population size as a function of time, averaged over 10^ simulation runs, for r values 
near the transition points at g = and g = 0.05 on the spatial lattice, with v = 1. Dashed and solid 
lines indicate the population growth for systems at dynamic percolation and directed percolation 
transitions respectively, showing that these transitions have the characteristic properties of those 
universality classes. (B) Critical r for the host-consumer ecosystem with v = 1. The vertical 
axis, where g = 0, maps onto the SIR epidemic process (described in [30]), while the line g = 1 
closely resembles the SIS dynamics seen in the epidemiological lattice model of [30] when rates for 
initial and secondary infections are the same. The transition line crosses over from the dynamic 
percolation universality class at 5 = to directed percolation between g = 0.015 and g = 0.02. 
Red Xs indicate the transition curve for the host-consumer dynamics played on a random regular 
graph (RRG) of uniform degree 4; the dashed line connecting them is to guide the eye. The RRG 
transition is neither directed percolation nor dynamic percolation. 

One can calculate the minimum sustainable value of r in the mean-field approximation 
[46] by balancing the birth and death rates. If the host population is small compared 
to the total ecosystem size, then the minimum sustainable r is the value which satisfies 
kr = V, where k is the number of neighbors adjacent to a site. For the parameters used 
in Figure 2(A), this value would be 0.05, a factor of 4 lower than that seen in the spatial 
evolving ecosystem. Figure 4(B) plots the minimum sustainable r in mean-field and spatial 
systems as a function of g. Both the minimum r seen in the scenario with mutation and 
the minimum sustainable r in the case without mutation are significantly larger than the 
mean-field threshold value. This discrepancy between the spatial and non-spatial scenarios 
further indicates the importance of incorporating spatial extent into epidemiological models. 

It is instructive to compare the minimum sustainable transmissibilities shown in Fig- 
ure 4(B), with the results of earlier investigations of spatial epidemiological models. For 
simplicity, we consider an SIR model in the case of rapid recovery, so that individuals which 
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are I[nfected] will be R[ecovered] in one simulation timestep. The infection cannot spread 
unless the disease transmissibility is sufficiently high [30, 66]. The minimum sustainable 
transmissibility is 0.5, the same as in the host-consumer model when the host growth rate 
g = 0. See Figure 4(B). Furthermore, the transition in the SIR lattice model is also in the 
dynamic percolation universality class, just as is the transition in the host-consumer model. 
If the SIR model is modified to allow reinfection, becoming an "SIRI" model, the mini- 
mum sustainable transmissibility decreases as g increases. The transition is, however, still 
a dynamic percolation one, up until the initial infection {S — )■ /) and reinfection {R — > /) 
probabilities are equal, in which case the transition becomes a directed percolation one in- 
stead. Contrast this with the host-consumer ecosystem, in which increasing the host growth 
rate g up from zero also decreases the critical value of the transmissibility — allowing the con- 
sumer population to sustain itself at lower r — but the transition is a directed percolation, 
not dynamic percolation, phenomenon. 

When the hosts reproduce to fill empty sites as quickly as possible {g = 1), then sites left 
open by dying consumers rapidly become available for consumer recolonization, and empty 
space does not play as significant a role. We expect, therefore, that ioi g = 1 the host- 
consumer model closely resembles the "SIS" epidemiological model, in which all non-infected 
individuals are intrinsically equally susceptible to infection [12, 30]. Indeed, we find that the 
transition points are the same for the two models: both have a directed percolation-class 
transition at a transmissibility of ^0.295 (see Figure 4(A) and reference [30]). 

Thus, in the limiting cases of g = and g = 1, the host-consumer model is equivalent to 
the SIR and SIS epidemiological models, while over the range in between, they differ. The 
key difference is that reoccupying an empty site with a consumer requires prior recolonization 
by a host, whereas the vulnerability of a R[ecovered] individual to becoming I[nfected] is 
defined as an intrinsic property of the R[ecovered] type. This changes the role of ecology: 
both models incorporate space, but the effect of spatial extent is different. 

C. Upper Threshold: Malthusian Catastrophe 

The upper cutoff in Figure 3 is the value of r above which the consumer population is 
unsustainable. Consider the limiting scenario where ~ 0. If hosts do not reproduce into 
available empty sites, we have an epidemic process. In this case, the minimum r is also a 
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maximum. Below r = 0.5, a consumer injected into a lattice of hosts will produce a consumer 
strain (which we can think of as an infection) which survives for a finite number of generations 
and then dies out, leaving the lattice filled with hosts (susceptibles) marred by a small patch 
of empty sites (recovered individuals). Above the transition point, a single consumer gives 
rise to an expanding wave of consumers which propagates over the lattice, leaving empty 
sites in its wake, until it consumes all the hosts in the ecosystem. This constitutes a process 
of annular growth [30]. No finite ecosystem can sustain annular growth indefinitely. If 
the host growth rate g is made non-zero, then hosts can recolonize sites left empty by the 
expanding consumer population, opening the possibility of host-consumer coexistence in an 
ecosystem of dynamically formed and re-formed patches. Figure 1 illustrates an example of 
this phenomenon. 

We can, therefore, interpret the upper cutoff on consumer sustainability as a Malthusian 
catastrophe due ultimately to the limited amount of available hosts [65]. For the case of a 
single consumer injected into a lattice otherwise filled with hosts, the cutoff is due to the 
limited size of the system, a type of outcome often called a finite-size effect. This is the key 
to understanding what happens when multiple types of consumer are present on the same 
lattice, and in particular the case we study in the next section, where an invasive consumer 
variety is introduced to an ecosystem where native hosts and consumers have already formed 
a dynamic patch distribution. The environment experienced by the invasive variety is that 
formed by the native species, and the "finite size" of the resources available to the invasive 
variety is not the size of the whole lattice, but that of a local patch [65]. 

III. TIMESCALE DEPENDENCE OF INVASION SUCCESS 

A key question about an ecological system is whether a new variety of organism, having 
a different genetic character and phenotypic trait values, can successfully invade a native 
population. One relevant measure is the expected relative growth rate of the number of 
offspring of a mutant individual within a native population, the relative growth rate of 
a mutant strain. This rate, known as the invasion fitness, is often used to investigate 
the stability of an evolutionary ecosystem [11, 67, 68]. If the invasion fitness is found to be 
positive, the native variety is judged to be vulnerable to invasion by the mutant. Conversely, 
if the invasion fitness is found to be negative, the native variety is deemed to be stable. We 
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shall see that for this host-consumer ecosystem this method gives qualitatively incorrect 
predictions for evolutionary dynamics. 

If a mutant consumer strain with fixed transmissibility successfully invade a 

population of transmissibility tq < Tm, then we expect the time-averaged value of r seen in 
the evolving system to be larger than tq. To investigate this, we simulate r = 0.45 mutants 
entering an ecosystem whose native population has r = 0.33. Initially, the mutants prosper, 
but they ultimately fail to invade. As shown in Figure 5, the probability of a r = 0.45 
strain surviving for tens of generations after injection is larger than that of a r = 0.33 
strain. That is, mutants with the higher r can out-compete the neutral case. However, after 
~ 74 generations, the survival-probability curves cross. Observed over longer timescales, 
the mutant strain is less successful than the native variety. This key result manifests the 
distinctive properties of the spatial structure of the model. The underlying reason for this 
result is that the mutants encounter the resource limitations imposed by the patchy native 
population. Over short timescales, the mutant strain enjoys the resources available within 
the local patch, consuming those resources more rapidly than can be sustained once it 
encounters the limitations of the local patch size. Short-term prosperity is not a guarantee 
of medium- or long-term success. 
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FIG. 5: (A) Survival probability as a function of time (logarithmic axes) for three scenarios: 
injecting mutants with the same transmissibility as the native consumers, injecting mutants with 
transmissibility r = 0.45 which is higher than the r of the native consumers, and injecting r = 0.45 
mutants with the swapping procedure described in §IV enabled. (The "background" was defined 
by 5 = 0.1, V = 0.2, r = 0.33.) Adding the swapping effect allows voracious mutants to persist 
long enough that they have a non-zero chance of taking over the ecosystem. (B) Detail of the plot 
showing the crossover point (linear axes). 
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This is to be contrasted with what happens in a well-mixed ecosystem. In the well- 
mixed scenario, consumer strains with higher r successfully invade and displace the native 
population with a high probability. The invasion success is consistent with the dynamics of 
a continuously evolving ecosystem. If r is made an evolvable trait in panmictic systems, the 
average r of the population will tend to 1. There is no difference in a well- mixed scenario 
between short-term and long-term success. This follows from the lack of distinction between 
local patches and large-scale structure. 

Random regular graphs (RRGs) have been used as approximations to incorporate the 
effects of spatial extent into population models, as they make for more tractable mathemat- 
ical treatments (including the "moment closure" methods discussed in the next section), 
although they are typically less realistic than spatial lattices [52]. In an RRG, each node 
has the same number of neighbors, as they do in a lattice network, but the connections are 
otherwise random. The network structure is set at the beginning of a simulation and does 
not change over time. When we simulate our host-consumer ecosystem on an RRG, we 
find that an invasive consumer strain with higher transmissibility r can out-compete and 
overwhelm a native consumer population with lower r. In one simulation run, using the 
native and invasive r values of 0.33 and 0.45 respectively, 2,233 out of 10,000 invasions were 
successful, whereas on the lattice no invasion succeeded using the same parameters. Thus, 
the RRG does not capture the essential features of the spatial scenario. In particular, our 
results show that the RRG case is more like the well-mixed case than the spatial lattice, as 
far as stability against invasion is concerned. 

We have found that spatial models evolve r values which are smaller than the maximum 
possible, as would be found in panmictic ecosystems. This means that consumer populations 
do not reproduce at the maximum rate. Whether predators can be "prudent" in their 
resource consumption has been a contentious issue in evolutionary biology [10]. Because the 
average r does not tend to 1, it appears that consumers in the spatial ecosystem exhibit 
a kind of "prudence". Yet, at any time, an individual consumer is as likely to spawn an 
offspring with higher r as it is one with lower r. (Barring, of course, the boundary cases 
r = and r = 1, beyond which the transmissibility cannot go.) Individual consumers do not 
judge the host resources available in their vicinity in order to adjust their reproduction rate 
accordingly. Instead, the apparent "prudence" arises from natural selection over the history 
of the population, eliminating unsustainable strains through local Malthusian catastrophies 
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many generations after they arise. The distribution of r across the population is, at any 
moment, the product of effects at work during the past, some of which required many 
generations to play out. 

Our investigation builds on earlier work which studied the timescale dependence of fitness 
indicators in spatial host-consumer ecosystems [44, 45]. In this paper we have augmented 
the prior work by considering the survival probability to show the effects of varying r. 
In addition, we reported the case of a mutant strain invading a background population, 
clarifying the conceptual and quantitative results of those earlier works, which considered 
instead scenarios complicated by multiple ongoing mutations. 

IV. PAIR APPROXIMATIONS 

The inadequacy of mean-field treatments of spatial systems motivates the development 
of more elaborate mathematical methods. In this section, we review one such methodology, 
based on augmenting mean-field approximations with successively higher-order correlations, 
and we test its applicability to our host-consumer spatial model. The numerical variables 
used in this methodology are probabilities which encode the state of the ecosystem and 
can change over time. One such variable might be, for example, the probability Pa that an 
individual chosen at random is an organism of type a. Another might be Pab, the probability 
that a randomly- chosen pair of neighboring individuals will have one member of type a 
and the other of type b. The change of these quantities over time is usually described 
by differential equations, for which analysis tools from nonlinear dynamics are available 
[11, 32, 65, 68-71]. 

The importance of the joint probabilities Pab is that they reflect correlations which mean- 
field approximations neglect. To understand the relevance of the joint probabilities Pab, note 
that if an invasive mutant variety forms a spatial cluster near its point of entry, then the 
average density of invasive mutants in the ecosystem, pm, will be low, while the conditional 
probability that a neighbor of an invasive individual will also be of the invasive type, qM\M = 
Pmm/vmi will be significantly higher. A discrepancy between qa\b and Pa can persist when 
the ecosystem has settled into a quasi-steady-state behavior, and is then an indicator of 
spatial pattern formation. Figure 6(A) shows an example, plotting qc\c measured in a host- 
consumer ecosystem against pc- The salient point is that in a well-mixed scenario (where 
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we expect the mean-field approximation to be applicable), the average consumer density pc 
and the consumer-consumer pairwise correlation qc\c are essentially equal over time. In the 
spatial lattice scenario, pc and qc\c are noticeably different. 

Treating the correlations qa\b as not wholly determined by the probabilities pa is a way 
of allowing spatial heterogeneity to enter an analytical model. Whether it is a sufficient 
extension in any particular circumstance is not, a priori, obvious. Typically, the differential 
equations for the pair probabilities Pab depend on triplet probabilities Pabc, which depend 
upon quadruplet probabilities and so forth. To make progress, we truncate this hierarchy 
at some level, a technique known as moment closure [58, 67, 69, 72]. Moment closures 
constitute a series of approximations of increasing intricacy [51, 73]. The simplest moment 
closure is the mean field approximation; going beyond the mean field to include second- 
order correlations but neglecting correlations of third and higher order constitutes a pair 
approximation. These approximations destroy information about spatial structure which 
may be necessary to account for real- world ecological effects [67]. 

Pair approximations have been used to test for the existence of an Evolutionary Stable 
Strategy (ESS) in a system — that is, a strategy which, when established, cannot be suc- 
cessfully replaced by another [74]. In addition to the limitations of pair approximation for 
representing patch structure [68] , as we saw in the previous section, the question of whether a 
mutant strain can initially grow is distinct from the question of whether that strain achieves 
fixation or goes extinct [37, 43, 75-79]. The former is a question about short-term behavior, 
and the latter concerns effects apparent at longer timescales. This distinction is often lost 
or obscured in analytical treatments. The reason is that one typically tests whether a new 
type can invade by linearizing the corresponding differential equations at a point where its 
density is negligible. However, this only reveals the initial growth rate (see the fixed-point 
eigenvalue analysis in [11, 67, 68]). 

In our model, we can directly test the efficacy pair approximations in a completely general 
way. The key idea is to transplant individuals in such a way that the variables used in the 
moment-closure analytical treatment remain unchanged. At each timestep, we look through 
the ecosystem for isolated consumers, that is, for individual consumers surrounded only by 
hosts and empty sites. We can exchange these individuals without affecting the pairwise 
correlations. For example, if we find a native-type consumer adjacent to three hosts and one 
empty site, we can swap it with an invasive- type consumer also adjacent to three hosts and 



17 



A 04 




B 0.16 




0.2 0.4 0.6 0.8 
Swapping probability 



1.0 



FIG. 6: (A) Pairwise conditional probability qc\c plotted against the average density of consumers, 
PC, for three variations on the host-consumer model: a well- mixed case in which mean- field theory 
is applicable, a random regular graph (in which each site has exactly four neighbors) and a 2D 
square lattice. The dotted line, pc = qc\c^ indicates the mean- field approximation. 10^ timesteps 
were computed for each case. The well-mixed case is simulated by dynamically rewiring sites at each 
time step, precluding the generation of spatial heterogeneity; consequently, the pairwise correlation 
qc\c is essentially pc {R^ = 0.953). The random regular graph (RRG) with random but static 
connections does develop spatial heterogeneity so that qc\c is not the same as pc = 0.581). 
However, the discrepancy is significantly stronger in the lattice case {B? = 0.304). (B) Success 
rate of invasive mutant strains as a function of swapping probability. Voracious mutant strains 
with T = 0.45 are introduced into a lattice ecosystem defined by a host growth rate of 5 = 0.1, 
a consumer death rate v = 0.2 (the same for both consumer varieties), and a native consumer 
transmissibility of r = 0.33. Average success rates are found by simulating 2000 invasions per value 
of the swapping probability parameter; error bars indicate 95%-confidence intervals. Increasing the 
fraction of possible swaps which are actually performed makes the voracious invasive strain more 
likely to take over the ecosystem. 



one empty site. We can also exchange isolated pairs of consumers in the same way. The 
variables used in the moment-closure treatment remain the same. Were the moment- closure 
treatment valid, we would expect to see the dynamics remain unchanged. 

When we perform the simulation, however, swapping does affect the dynamics, and quite 
strongly. With this type of swapping in effect, mutants with r = 0.45 can invade a native 
population with r = 0.33, doing so 1,425 times out of 10,000 injections in one simulation 
run. Without swapping, the number of successful invasions is zero. 

Swapping creates a new ecosystem model whose moment-closure treatment is the same as 
of the original. The behavior of invasive strains is different, because transplanting organisms 
allows invasive varieties to evade localized Malthusian catastrophies. Swapping opens the 
ecosystem up to invasive strains, since, in essence, it removes individuals from the "scene of 
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the crimes" committed by their ancestors. 

This type of swapping is, to our knowledge, a new test of moment-closure validity. Ran- 
domized exchanges have been incorporated into computational ecology simulations for differ- 
ent purposes. For example, research on dispersal rates in an island model shuffled individuals 
in such a way that the population size of each island was held constant [80]. 

If, instead of performing every permissible swap, we transplant organisms with some 
probability between and 1, we can interpolate between the limit of no swapping, where 
invasions always fail, and the case where pair approximation is most applicable and invasions 
succeed significantly often. The results are shown in Figure 6(B) and indicate that the 
impact of swapping becomes detectable at a probability of !^0.25 and effectively saturates 
at a probability of ^0.9. 

It is instructive to compare the spatial lattice ecosystem with the host-consumer model 
defined on a random regular graph (RRG). Consider, first, the pairwise correlation qc\c- On 
an RRG, the underlying network topology provides enough of a sense of locality that pc 
and qc\c are less strongly correlated than in the well- mixed case. See Figure 6(A). We can 
also implement swapping on the RRG, where invasions can succeed without it; as expected, 
swapping does not affect the success rate on the RRG. (With 10,000 simulated invasions 
for each case, the 95%-confidence interval for the difference in success rates between full 
swapping and none is 0.004 ±0.01.) 

Other numerical investigations have also suggested the limited utility of pair approxi- 
mations [52, 65, 69, 70, 74, 81, 82]. One reason cited is the problem of phase transitions. 
Moment closures fail when a system is near a transition such as the directed percolation 
and dynamic percolation ones described above [51]. However, the timescale dependence of 
reproductive fitness, which stems from non-equilibrium pattern formation, indicates that 
these problems are generic, occurring even away from the phase transition. 

Our analysis implies that pair approximations are inadequate for analysis of systems 
with spatial inhomogeneity. One might ask whether including triple and other higher-order 
corrections could help. However, it is known that this series approximation is non-analytic 
at phase transitions [51]. Non-equilibrium pattern formation will necessarily also be poorly 
described, at least until the order of expansion reaches the characteristic number of elements 
in a patch. Given the algebraic intricacy of higher-order corrections to pair approximations 
[52, 65, 81], it is useful to know in advance whether such elaborations have a chance at 
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success. If approximation techniques based on successively refining mean-field treatments 
are blind to important phenomena, then we should build our analytical work on a different 
conceptual foundation. 

V. CONCLUSIONS 

Understanding the effects of spatial extent is a vital part of evolutionary ecology. Spatial 
extent changes the quantitative and qualitative characteristics of a model's evolutionary 
behavior, compared to well-mixed models. The short-term success rate of novel genetic 
varieties is not indicative of their long-term chance of success relative to the prevalent type. 
Standard stability criteria fail to reflect the actual stability achieved over time. We must 
instead consider extended timescales because they are determined by spatial patterns, whose 
ongoing formation is an intrinsic part of nonequilibrium evolutionary dynamics. Our analysis 
shows that there are dramatic differences between spatial models and mean-field models, 
which simplify away heterogenity through mixing populations, averaging over variations or 
mandating a globally connected patch structure. We have further shown that transplanting 
organisms dramatically changes the dynamics of spatial systems, even when we preserve 
local correlations as would be considered in a pair approximation treatment. Our results 
prove that any model striving to capture the effects of heterogenity that does not change its 
behavior with organism transplanting cannot fully capture the dynamics of spatial evolution. 

In our host-consumer models, each individual either survives or it does not, and any 
individual can produce zero or more offspring; that is to say, an "individual fitness" (in 
the terminology of [83]) is a well-defined concept. To find expected individual fitness, or 
average individual fitness, we must define a set of individual organisms over which to take 
an average, which is the very concept we have established to be problematic. Consequently, 
derived notions of fitness, which depend on comparisons between such averages [83] , become 
elusive, context-dependent quantities. This dependence is both temporal and spatial: Av- 
erage relative fitness in one generation is not necessarily a good measure of the long-term 
success of a strain in one, or a combination of, the broad variety of dynamically-generated 
niches. This problem is not the same as the traditional concept of variation of fitness across 
a static set of niches, because the niche dynamics ensures that evolutionary outcomes are 
not reflected in any reasonable definition of the average. 
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One might be inclined to call the behavior of voracious invasive consumers "frequency- 
dependent fitness" [4, 5], as the invasive strain is successful initally when rare but fails when 
it becomes more common. The term "frequency-dependent fitness" is, however, a misnomer 
in this context. The organism type is rare and successful when it is newly introduced, but 
as it declines to extinction it becomes rare and unsuccessful. Frequency, being defined by 
an average over the whole ecosystem, is only a proper variable to use for describing the 
ecosystem in the panmictic case. One might refine the concept of global frequency to the 
notion of local frequencies, which can vary from place to place across the ecosystem. Here, 
though, even the extension to local frequencies fails to capture the evolutionary dynamics. 
Including local frequencies implies incorporating local correlations or local environmental 
characteristics. However, the breakdown of moment-closure techniques implies that defining 
fitness as a function of organism type together with average local environment [84] will, in 
many circumstances, not go far enough. 

We find that we cannot assign a meaningful invasion fitness value to an invasive variety 
of organism. Consequently, we cannot assign a fitness value to a genetic characteristic which 
defines such an organism type. That is, defining a trait or allele fitness is likewise invalid. To 
restate the host-consumer model in genetic terms, we could say that an individual consumer 
can have one of two possible alleles of the "transmissibility gene", one coding for r = 0.33 
and the other for r = 0.45. A mean-field treatment would then involve specifying the 
fraction of the population which carries the r = 0.33 allele versus the fraction which carries 
the T = 0.45 variant. Such a description, however, is too oversimplified to let us predict what 
will happen in the ecosystem. All of the issues we identified imply that the evolutionary 
dynamics cannot be characterized using the allele frequencies at a particular time. 

If we can no longer summarize the genetic character of a population by an allele 
frequency — or a set of allele frequencies for well-defined local subpopulations — then com- 
puting the fitness of a genotype from its generation-to-generation change in frequency is a 
fruitless task. In a world which exhibits nonequilibrium spatial pattern formation, allele 
frequencies are the wrong handle for understanding the dynamics of natural selection. For- 
mally, the conventional assumption that the allele frequencies are a sufficient set of variables 
to describe evolutionary dynamics is incorrect. The spatial structure itself is a necessary 
part of the system description at a particular time in order to determine the subsequent 
generation outcomes, even in an average sense. 
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The timescale-dependence issues which arise in spatial host-consumer ecosystems exist in 
a wider context. Multiple examples indicate that initial success and eventual fixation are only 
two extremes of a continuum which must be understood in its entirety to grasp the stability 
of a system. In the study of genetic drift, it has been found that neutral mutations can 
fixate and beneficial mutations fail to fixate due to stochasticity [4]. Likewise, in the study 
of clonal interference [78] , one beneficial mutation can out-compete another and prevent its 
fixation. Furthermore, classical genetics makes much use of the Price Equation for studying 
the change in a population's genetic composition over time [6, 7], and it is well known that 
analytic models built using the Price Equation lack "dynamic sufficiency". That is, the 
equation requires more information about the current generation than it produces about the 
next [5, 6, 85], and so predictions for many-generation phenomena must be made carefully, 
if they can be made at all. Modeling approaches which are fundamentally grounded in the 
Price Equation, such as "neighbor-modulated" fitness calculations [7, 19, 41, 86] and their 
"multilevel" counterparts [7, 86, 87], are not likely to work well here, as the analyses in 
question consider only the short-timescale regime. In addition, those particular analyses 
which address host-consumer-like dynamics either rely on moment closures [41] or they 
assume a fixed, complete connection topology of local populations which are internally well- 
mixed [19, 87]. These simplified population structures are quite unlike the dynamical patch 
formation seen in the host-consumer lattice model. (Wild and Taylor [88] demonstrate an 
equivalence between stability criteria defined via immediate gains, or "reproductive fitness" , 
and criteria defined using fixation probability; however, their proofs are explicitly formulated 
for the case of a well-mixed population of constant size, neither assumption being applicable 
here.) 

Our results also have significance in the context of adaptive-network research. The field 
of adaptive networks studies scenarios in which the dynamics of a network and the dynamics 
on that network can occur at comparable timescales and feed back on one another. When 
building such a model, one chooses a rewiring rule which specifies how the states of the 
network nodes affect the changes of the link structure which connects them. As discussed 
earlier, in a spatial host-consumer system, connections can form between resource patches 
when hosts reproduce into intervening empty sites, whereas in adaptive-network studies of 
SIR models, new connections are made without regard to spatial proximity. The difference 
we have seen between lattice behavior on one hand and RRG or swapping-enabled behavior 
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on the other emphasizes the need to study the effect of spatial proximity on hnk rewiring. 
While the structure-erasing nature of unconstrained rewiring among susceptible hosts has 
been acknowledged [61, 62], new rewiring rules which reflect spatial and community structure 
have yet to be systematically investigated. "Myopic" rewiring rules, such as restricting the 
set of eligible new partners to the neighbors of a node's current partners, have on occasion 
been considered, but in contexts other than epidemiology, like evolutionary game theory 
[89, 90], making the endeavour of exploring such rules in epidemic models all the more 
worth pursuing. 

Our analysis of transplanting organisms can be considered parallel to real world concerns 
and manifest effects of invasive species [37, 42]. These are among the well-established ex- 
amples of situations in which spatial extent influences evolutionary dynamics [79, 91-94]. 
Identifying specific implications of the issues explored in this paper for particular biological 
systems [39, 43, 48, 49, 79, 93, 95] requires field and laboratory work. 
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